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Abstract. Multimodal oncological strategies which combine chemother- 
apy or radiotherapy with hyperthermia have a potential of improving the 
efficacy of the non-surgical methods of cancer treatment. Hyperthermia 
engages the heat-shock response mechanism (HSR), which main compo- 
nent (the heat-shock proteins) is known to directly prevent the intended 
apoptosis of cancer cells. Moreover, cancer cells can have an already par- 
tially activated HSR, thereby hyperthermia may be more toxic to them 
relative to normal cells. On the other hand, HSR triggers thermotoler- 
ance, i.e. the hyperthermia treated cells show an impairment in their 
susceptibility to a subsequent heat-induced stress. This poses a ques- 
tion about the efficacy and about an optimal strategy of the therapy 
combined with hyperthermia treatment. 

We adapt our previous HSR model and propose its stochastic extension, 
which we then analyse using the approximate probabilistic model check- 
ing (APMCl technique. We formalise the notion of the thermotolerance 
and estimate the intensity and the duration of the HSR-induced thermo- 
tolerance. Finally, we quantify the effect of a multimodal therapy based 
on hyperthermia and a cytotoxic effect of bortezomib, a proteasome in- 
hibitor, and we propose an optimal strategy for combining these two 
modalities. 

By mechanistic modelling of HSR we are able to support the common 
belief that the combination of cancer treatment strategies increases ther- 
apy efficacy. Moreover, our results demonstrate feasibility and practical 
potential of | APMCl m analysis of stochastic models of signalling path- 
ways. 



1 Introduction 

Most of the non-surgical methods of cancer treatment (e.g. chemother- 
apy and radiotherapy) are based on the principle of putting some kind of 
stress on cancer cells to induce their death. Unfortunately, in many cases 
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the above methods fail. The fact that heat-shock proteins (HSP) prevent 
apoptosis induced by different modalities of cancer treatment explains 
how these proteins could limit the application of such anti-cancer ther- 
apies. In order to improve the efficacy of these treatments, some effort 
is focused on the multimodal oncological strategies which usually com- 
bine treatment of chemotherapy and hyperthermia, or radiotherapy and 
hyperthermia. 

1.1 Heat-shock response in cancer treatment 



HSP are a group of highly conserved proteins involved in many phys- 
iological and pathological cellular processes. They are so called chaper- 
ones, as which they help with folding of new and distorted proteins into 



their proper shape. In principle, HSP synthesis increases under stress 



conditions. Subsequently, induction of |HSP| increases cell survival and 



stress-tolerance. Elevated expression of different members of HSP family 
in tumour cells has been detected in several cases. Despite its importance, 
little is still known about how exactly [HSP are involved in different pro- 
cesses related to cancer development. 

In this work we're interested in the heat-shock inducible isoform of 
heat-shock proteins 70 kDa [Hsp70; see, e.g., [27]. For a sake of clarity, 
we will denote Hsp70 protein by HSP, and use the former only if context 
might be unclear. 

Hyperthermia is a therapeutic procedure used to raise the tempera- 
ture of a region of the body affected by cancer [29] . Body tissue is, globally 
or locally, exposed to temperatures up to 45° C. Moderate hyperthermia 
treatment maintains temperatures for about an hour in a 40-42°C range, 
in which most of normal tissues are not damaged. Besides characteris- 
tics specific to cell type, the effectiveness of hyperthermia depends on the 
temperature achieved during the treatment, as well as on the length of the 
treatment [29| 9] . Currently, hyperthermia effectiveness is under study in 
clinical trials, including combination with other cancer therapies [29, 30J. 

A synergistic interaction of radiotherapy and hyperthermia as well 
as some cytotoxic drugs and hyperthermia has already been confirmed 
in experimental studies [9]. In particular, Neznanov et al. [T7] demon- 



strated, in vitro, that induction of HSR by hyperthermia enhances the 



efficacy of a proteasome inhibitor called bortezomib — a U.S. -approved 
drug for treatment of a multiple myeloma and mantle cell lymphoma 



16] . Basically, hyperthermia engages the HSR mechanism which main 



component are the anti-apoptotic HSP Cancer cells have an already par- 



tially activated HSR because they are constitutively coping with higher 
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level of misfolded protein (mainly because of their rapid rate of prolifer- 
ation as well as specific intracellular conditions). Therefore, in principle, 
sufficiently increased level of misfolded proteins, as obtained by, e.g., hy- 
perthermia, can not be matched by cell's HSR capacity and, in effect, 
such enhanced proteotoxic stress can be more toxic to them relative to 
normal cells |17j . 

On the other hand, after a heat-shock, all cell types show an im- 
pairment in their susceptibility to heat-induced cytotoxicity. This phe- 
nomenon, known as thermotolerance, is triggered by HSR [it is at least 
partially based on the induction of HSP; 9]. Thermotolerance is, in princi- 
ple, reversible and persists for usually between 24 and 48 hours [29] . Due 
to this phenomenon the applicability of the combined hyperthermia ther- 
apy may be, counter-intuitively, initially limited. This poses a question 
about the actual efficacy and about an optimal strategy of the hyperther- 
mia treatment. 



1.2 Related HSR models 



There are basically two types of approaches to HSR modelling. They 
differ in how the heat-induced protein misfolding is captured in the model. 
In the first approach, the temperature is implicitly represented by a model 
variable and its increase is modelled by a direct change of the system 
state. More specifically, when temperature change is modelled by addi- 
tion of an intermediate species, which causes native proteins to misfold 
[see, e.g., [22]. In the second approach, temperature is explicitly expressed 
in model parameters. To that end, Peper et al. |20j proposed an integral 
model of HSP regulation, where temperature-dependence is expressed in 
kinetic rate constant of native proteins misfolding reaction; the temper- 
ature change is modelled by a direct change of kinetic parameter of one 
reaction. This approach, with the same temperature-dependence function 



has been followed up in some of the recent papers which include HSR 
modelling (26] [2TJ [15] • In between is an approach presented in Rieger 
et al. [23] is based on a model reduction by excluding the intermediate 
species (a kinase of misfolding) , thus, conceptually being a case of the first 
approach, but formally being a case of the second approach. More specif- 
ically, Rieger et al. [23] modelled temperature perturbation by a change 
of kinetic rate constants. However, in his approach misfolding is based on 
enzyme kinetics approximation of Michaelis and Menten |13| with arbi- 
trarily selected kinetic values for four different temperatures]. Contrary 
to that, Peper et al. J2U| modelling approach follows purely mass action 
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kinetics and the misfolding rate is based on experimental data on the 
level of protein denaturation for a range of 35-50°C [12, Fig. 8B]. 

As far as the stochastic HSR models are considered, the counterpart 
of Petre et al. [21] model was presented by Mizera and Gambin |15) . 
The authors thoroughly validated deterministic approximation. Another 
stochastic and deterministic pair of HSR models was recently presented 
in Proctor and Lorimer |22j . 

For our case study we have chosen the explicit temperature-dependent 
denaturation rate approach (second one), in which the model proposed 
by Szymahska and Zylicz [26] is the simplest. With many simplifications, 
such as direct trimerisation or no misfolding of |HSR components, this 



model is able to capture the overall qualitative behaviour of the HSR 



mechanism. We would like to emphasise that, according to authors knowl- 
edge, none of the above-mentioned models have been used to investigate 
either thermotolerance or multimodal oncological strategies. 



1.3 Probabilistic model checking of biochemical pathways 



A classical modelling paradigm of biological systems, biochemical path- 
ways in particular, uses differential equations to define the evolution of 
average molecular concentrations over time. In this paper we adopt a 
stochastic approach, where the evolution of individual molecules is mod- 
elled, with rates of interactions controlled by exponential distributions — 



a continuous-time Markov process (CTMC). This approach allows for 



computing exact quantitative measures as opposed to average values ex- 
pressible in the classical approach. 



The probabilistic model checking (PMC) is an extension of model 



checking, a well-established formal verification technique successfully ap- 
plied in many of analyses of computer systems. Adopting a stochastic 
modelling paradigm allows us to take advantage of a efficient probabilis- 
tic model checker PRISM. The motivation for using PMC] is the belief 
that when used in conjunction with other, well-established approaches, 
such as ODE numerical solutions, PMC may offer a greater insight into 



the complex interactions present in biological pathways. 

Recent research demonstrated a considerable success in adapting [PMC] 
to analysis of signalling pathways models. Among such systems analysed 
recently with the use of PRISM are: the MAPK cascade [II], the RKIP 
inhibited ERK pathway [3j, the FGF pathway [8] or the T cell signalling 
pathway [TH]; Calder et al. [3] contains an overview. 
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1.4 Our results 

The main purpose of this work is to contribute to the understanding 
of the involvement of the IHSRI mechanism in multimodal cancer ther- 
apies. To this end we use a refined version of the deterministic model 
of Szymahska and Zylicz |26j . which despite of its simplicity provides a 
sound qualitative model of the HSR] mechanism. In this work, we develop 



a stochastic version of this model, represented by chemical master equa- 



tion (CME) or, equivalently, CTMC, which we then analyse using the 



PMC technique. To ensure the feasibility of this approach we have used 



approximate PMC] techniques, implemented recently in the PRISM tool 



[see, e.g., HE]. 

Next, we formalise the notion of thermotolerance i.e. the memory of 
the system of the previous temperature perturbation, or system desensi- 
tisation with respect to the second consecutive heat-shock. We compute 
the intensity and the duration of the |HSRf induced thermotolerance. Fi- 
nally, we give a proof of concept quantification of an effect of a combined 
therapy of hyperthermia and an bortezomib-induced proteasome inhibi- 
tion. Based on that, we propose an optimal strategy for combination of 
heat-shock and the inhibitor. In principle, our results support the com- 
mon belief that the combination of cancer treatment strategies increases 
therapy efficacy 



2 Model 



We consider the dynamics of synthesis of HSP and its interactions 
with key intracellular components of [HSR i.e.: HSP| the heat-shock fac- 
tor (HSF) and its trimer, which is a HSP transcription factor; HSP 
substrate — mainly denatured, misfolded native proteins; |HSP| gene - 



heat-shock element (HSE); and HSP mRNA. Figure [T] depicts the model 
scheme, Table[T]gives reaction list, whereas Table|2]gives the implied mass 
conservation constraints. Structural changes with respect to the previous 
version of this model by Szymahska and Zylicz [26] include explicit native 



protein species (Eqs (r9) and (rlO)) and separate HSP mRNA degradation 



(Eqs ( [rTT] ) and ([rig}). 

Tables [3] and [4] summarise, respectively, the variables of the model 
with their initial values, as well as description and values of all kinetic 
parameters. Parameters kj, and lj denote the mass action reaction rates 
constants, respectively, for the forward and backward reactions. Value of 
the reaction rate constant k^ depends on the given temperature T; more 
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Fig. 1: Scheme of the HSR model. On the left side of the scheme, the denaturation of native 
proteins P and denatured proteins S (substrate) refolding or degradation moderated by the 
|HSP| chaperones. On the right side, the adaptive |HSP| production loop, stimulated by |HSF| 
which trimerise and init iate |HSE transcription and HSP mRNA translation (dotted arrow). 
As a negative feed back, |HSP| molecules promote |HSF| t rimers dissociation and inhibit single 
|HSF| molecules by direct binding. To close the loop, in hi bited |HSF| molecules are forced out 
of the complex with HSP by the inflowing substrate. 



specifically: 

k£ « m u X (l - J£i_) X 1A AT mm" 1 , (1) 

where AT = T — 37 is a heat-shock temperature delta and mu = 0.0105 is 
a time-scale dependent multiplier. This function was originally proposed 
by Peper et al. [2D] , and recently reused in several other works [26] |2"T||15|. 
Eq. ([TJ is valid in a local temperature range, i.e. for T £ [37,45] °C |20j . 



2.1 Mathematical models 

Deterministic framework 

As mentioned, all reaction channels follow the mass action kinetics. 
The i th substrate amount in the rate of j th reaction is raised to the power 
of its' stoichiometric coefficient riij, where N = (n^) £ N 10 x N 16 denotes 
stoichiometry matrix of total of ten species involved in sixteen reactions, 



as described by Eqs (rl)-(rl2). 

In context of the deterministic modelling framework, state of our sys- 
tem is represented by the time dependent state vector y(i) of concentra- 
tions of reacting species yi (i = 1 . . . 10). Given the initial state yo (see 
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HSP:HSF + S 

3- HSF 
HSF 3 + HSE 

HSE:HSF 3 
HSP + HSF 3 
HSP + S 

HSP + HSF 

HSP 
HSP:S 

P 

mRNA 
mRNA 



k 3 
k 8. 



12 

k 9. 



HSP:S + HSF, 

HSF 3 , 
HSE:HSF 3 , 

HSE:HSF 3 + mRNA, 
HSP:HSF + 2 • HSF, 
HSP:S, 

HSP:HSF, 

0, 

HSP + P, 

S, 

mRNA + HSP, 



(rl) 

(r2) 
(r3) 

(r4) 
(r5) 
(r6) 

(r7) 

(r8) 
(r9) 

(rlO) 
(rll) 
(rl2) 



Table 1: The HSR 



biochemical reactions network. There are 12 reactions flrlh— ( rl2 1 , 4 



of which are reversible, making it 16 reactions in total; li (2 = 1,2,6,7) denotes reverse 
reaction rate constant. The T superscript denotes a temperature dependence. 



Ptot 


= P(t)+Stot(t), 






where S to t(t) = S(t) + HSP:S(t), 


(cl) 


HSF to t 


= HSF(t) + HSP:HSF(t) + 3 ■ HSF 3 (t)+ 






+ 3-HSE:HSF 3 (t), 


(c2) 


HSEtot 


= HSE:HSF 3 (t) +HSE(t). 


(c3) 



Table 2: Mass conservation laws in the iHSRl biochemical reactions network. Neither con- 
centrations nor molecules numbers are annotated because Eqs (|c"Tj)-(|c3"l) hold for all t > 
in both deterministic and stochastic models. 
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Description 


Value 


HSP 


free HSP 


0.309 


HSF 


HSF monomers 


0.151 


S 


substrate (denatured/misfolded protein) 


0.113 


HSP:HSF 


HSP:HSF interacting complexes 


2.588 


HSP:S 


HSP:substrate complexes 


1.126 


HSF 3 


active HSF (trimer form) 


0.044 


HSE 


free HSE 


0.957 


HSE:HSF 3 


bound HSE 


0.043 


mRNA 


HSP mRNA 


0.115 


P 


native proteins (denaturation susceptible) 8.761 



Table 3: Description of variables used in the HSR model, together with initial conditions 



for the deterministic version. These are a steady state concentrations in a non stressed cells 
(T = 37°C; see text for detail). Values are given in arbitrary scale molar concentration 
|[al]V[). 





Rate constant of. . . 


Value 


ki 


HSP:substrate association 


1.47 


li 


HSP:substrate dissociation 


0.0175 


k 2 


HSP:HSF association 


1.47 


la 


HSP:HSF dissociation 


0.0175 


k. 


HSF trimers association (activation) 


0.0805 


Is 


HSF trimers dissociation (inactivation) 


0.020125 


k 4 


HSP translation 


0.1225 


k 6 


HSP mRNA degradation 


0.0455 


k 6 


HSP:HSF dissociation and HSP:substrate association 


0.0805 


le 


HSP:substrate dissociation and HSP:HSF association 


0.00126 


k 7 


HSE:HSF 3 association 


0.1225 


It 


HSE:HSF 3 dissociation 


0.1225 


k 8 


HSP mRNA transcription 


0.1225 


k 9 


HSP degradation 


0.0455 


kio 


misfolded protein refolding (substrate degradation) 


0.049 


kS 


protein misfolding (substrate production) 


Eq. p) 



HSR model parameters. All values, except for k 5 are taken directly fro m th e Szy- 



Table 4: 

mahska and Zylicz [26J model and scaled for the experimental data fit (cf. Figure S.2 1. Note 
that all mass action rate constants (i.e. ki and lj) have units min -1 (|a.s.l V0~ ran < -" R ' +1 , 
where rank (R) is a rank of reaction R, i.e., a sum of left-hand side stoichiometric coeffi- 
cients. 
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Table [3]), the dynamics of the system is governed by a set of 10 ordinary 
differential equations (ODE), called reaction rate equations (RRE): 



dy(t) 

dt 



f(y(t)) = JVv(y(t)), 



(2) 



where v(y(t)) € M 16 is a vector of reaction rates at a time point t. Figure[2] 
depicts behaviour of this model, which starts in the state of homeostasis 
(i.e. in a steady state for T = 37°C), in response to the T = 42°C heat- 
shock. Amounts of species are arbitrarily scaled, for each of them sepa- 
rately, to obtain values of a similar order of magnitude (denoted a.s.M). 
Model is calibrated with respect to the HSE:HSF3 42°C experimental 



data (cf. Figure S.2). 



Stochastic extension 

In the stochastic approach we provide a set of linear, autonomous 



ODE one for each possible state of the system. Such set of equations 



is called the chemical master equation (CME). The solution of the fc-th 
equation at a time point t corresponds to the probability of the system 
being in that particular state at that time t. The system state is encoded 
as a vector X(i) G N 10 containing molecule numbers of all species. The 
j-th reaction changes the molecule numbers by: 

X(t) ->X(t) + nj, 

where nj is the j-th column of N. Denote by Oj(X(i)) the propensity 
function associated with the j-th reaction. Probability of j-th reaction 
taking place in the infinitesimally small time interval [t, t + dt) is accu- 



rately approximated by a,-(X(i))cft. CME describing the time-dependent 



distribution P(x, t) of system states, i.e. the probability that X(i) = x, 
is given by: 

= £ a . (x _ n . )P(x _ t) _ £ aj( x)P(x, t). (3) 

3 3 

Because the propensity functions a,- are time- independent, the CME rep- 



resents a continuous-time Markov process (CTMC). More precisely, the 
solution of the CME gives the transient probabilities for all states of the 
CTMC. 

For a stochastic model we used the scaling coefficient 5 which relates 
concentrations in the deterministic model to number of molecules in the 
stochastic model. Value of 5 corresponds to a number of molecules per 
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Fig. 2: Numerical simulations of the HSR||RRE model for a constant 42°C heating strat- 



egy. Simulation starts at a 37°C steady state. The upper plot depicts [HSP| response to 
the temperature-stimulated inflow of denatured proteins S (substrate). Free substrate is 
instantaneously bound into a HSP:S complex. Insufficient amount of free |HSP"1 causes its 
extraction from the HSP:HSF complex, forming an initiative response of the cell. Released 
in exchange |HSF| induces adaptive production of |HSP| molecules to complement its defi- 
ciency as indicated by accumulation of S, with peak at ca. 25 min. After over 120 min the 
excess of new HSP is used to inhibit HSF activity. System stabilises after ca. 450 min (not 
shown) with most of constantly inflowing S secured in the HSP:S complexes. The lower plot 
depicts the adaptive |HSP| production, stimulated by |HSF| |HSF|tr imerises and initiate [H5E| 
transcription to mRNA, followed by further translation to HSP| as visible by the shifted 
activity of subsequent components. 
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one unit of concentration, i.e. 5 ■ [S] = #5. This approach is equiva- 
lent to considering approximate stochastic model of packs of ■ \ V\/S 
molecules instead of single molecules (here is Avogadro constant and 
|V| is the solution volume). Given HeLa cell volume of circa 2500 /j,m s = 
2.5 • 10- 12 dm 3 [H], we have that for 5 ~ 1.5 • 10 12 the stochastic model 
is accurate. We adjust reactions constants accordingly, with -\V\ ■= 5 
[cf.[5]. 

Mean comparison 



We find that for 5 = 100 the RRE model is in a good agreement 



with the stochastic variant for both 37°C and 42°C. Visual comparison of 
ODE and stochastic simulations is presented in Figure [3| Table [5] presents 



comparison of stochastic mean value with RRE values for 5 equal to 100 
and 1000. Sources of stochastic model error are two-fold: the rounding 
errors due to the molecules packaging and the propensity constants ap- 
proximations, especially for the only reaction with rank (R) > 1, i.e. the 



HSF trimerisation (Eq. r2). Although for 5 = 1000 the relative error in 
steady state is ca. 10 times lower, the stochastic simulation paths, and 
thus their running time, are over 10 times longer (data not shown). We 
find 5 = 100 to be a good compromise between accuracy and efficiency 
for our proof-of-concept case study. 



Relative error +/- 95% ■ 



interval in 



% 



Species 


S-- 


Homeostasis 
= 100 5=1000 


5 


Heat-shock 
= 100 (5=1000 


HSP 


12.5 


+/- 


.68 


1.31 


+/- 


.19 


8.4 


+/- 


.55 


0.83 


+/- 


.16 


HSF 3 


12.1 


+/- 


.86 


1.45 


+/- 


.26 


9.4 


+/- 


.75 


0.71 


+/- 


.23 


HSP mRNA 


12.1 


+/- 


.79 


1.23 


+/- 


.24 


8.8 


+/- 


.67 


0.79 


+/- 


.21 


HSE:HSF 3 


11.4 


+/- 


.87 


1.37 


+/- 


.26 


8.5 


+/- 


.74 


0.86 


+/- 


.23 


HSF 


6.9 


+/- 


.72 


0.88 


+/- 


.24 


5.1 


+/- 


.68 


0.76 


+/- 


.23 


substrate 


2.5 


+/- 


.69 


0.34 


+/- 


.22 


2.9 


+/- 


.44 


0.33 


+/- 


.14 


HSP:HSF 


1.6 


+/- 


.06 


0.21 


+/- 


.02 


1.8 


+/- 


.09 


0.21 


+/- 


.03 


HSE 


0.6 


+/- 


.04 


0.06 


+/- 


.01 


0.5 


+/- 


.05 


0.05 


+/- 


.01 


HSP:substrate 


0.1 


+/- 


.17 


0.04 


+/- 


.05 


0.1 


+/- 


.05 


0.01 


+/- 


.02 



RRE 



Table 5: Estimates of a relative error of each species mean value with respect to its 
value, i.e. |E (#5) — [.S*] | / ; values are given in percent. Relative errors were calculated 
in homeostasis (T = 37°C) and the heat shock steady state (T = 42°C), for two scaling 
coefficient 5 values. Species are sorted according to error values in homeostasis for 8 — 100; 
from the least to the most c onsistent with the |RRE] solutions. Steady state mean values 



were estimated using APMC with 10 independent simulation samples for each species. 
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Fig. 3: Comparison of the stochastic simulations with respect to ODE numerical solutions 
for the HSR model. Both homeostasis (left column) and heat-shock (right column) con- 
ditions are compared. Each plots shows 10 sample stochastic trajectories, estima ted mean 
+/— standard deviation of a sample of 10 3 stochastic simulations, and an ODE numerical 
solution (black). Here, we assumed 100 molecules per unit of concentration, i.e. S = 100. 



Level of stochastic noise 

The variance-to-mean ratio: 

Var (X) 



VMR (X) 



E(X) 



quantifies noise of a species amount variable X = #S at a fixed time 
point in the stochastic model, with respect to the Poisson birth-death 
process [see, e.g., [2H]. Table [6] contains estimated steady state values of 
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VMR in our stochastic HSR model. These VMR values are significant 



for some of the crucial species, both for the state of homeostasis and the 
steady state during the heat-shock. 



VMR +/- 95% confidence interval 



Species 


Homeostasis 


Heat-shock 


HSP 
HSF 

HSP mRNA 


3.05 +/- 0.65 
2.41 +/" 0.35 
1.68 +/- 0.29 


3.14 +/" 0.74 

2.21 +/- 0.40 

1.60 +/- 0.34 


substrate 

HSE:HSF 3 

HSP:substrate 

HSF 3 

HSP:HSF 

HSE 


1.19 +/- 0.24 

0.81 +/- o.i2 
0.78 +/- 0.55 
0.78 +/- o.i2 

0.27 +/- 0.46 

0.10 +/- o.ii 


2.32 +/- 0.53 
0.85 +/- o.i4 
0.57 +/- 0.77 
0.79 +/- o.i5 
0.38 +/- o.60 
0.03 +/- o.i3 



Table 6: Estimates of VMR for each species in homeostasis (T = 37°C) and the heat-shock 
steady state (T = 42°C). VMR estimates were calculated for 5=100. Species are sorted 



according to the |VMR"1 values in homeostasis; from the most to the least disperse. Dashed, 
vertical line separates the over-dispersed and under-dispersed variables, with respect to the 
Poisson distribution. The dispersion doesn't change much with temperature, except for 



the substrate (highlighted). Mean and variance values were estimated using APMC with, 
respectively, 10 4 and 5 ■ 10 4 independent simulation samples for each species. 



The steady state amount of substrate, HSP, HSF and HSP mRNA is 
over-dispersed with respect to the Poisson distribution, indicating their 
high stochasticity in our model. In general noise of the species amounts in- 
creases for the higher temperature parameter value: mean VMR in home- 
ostasis is 1.23, whilst in the 42°C heat-shock 1.32 (ca. 7.5% higher; see 
Table [6| . This is only due to the almost two- fold increase in the substrate 
noise (highlighted) . Of course this may be attributed to the significant in- 
crease of the substrate amount during the heat-shock steady state. Note 
however, that, for instance, the number of HSP:S complexes is also sig- 
nificantly larger in the heat-shock steady state (cf. Figure [2]), and yet its 
noise is much smaller (ca. 27% less). 



3 Results 

3.1 Quantification of the thermotolerance phenomenon 

Thermotolerance can be described as a desensitisation with respect 
to a consecutive heat-shock, compared to the response to the first heat- 
shock. In other words, thermotolerance represents a memory of the system 
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about the first two, "on" and "off" temperature perturbations, leading to 
a decreased response to the subsequent "on" perturbation. Such system's 
memory is created by a propagating shift in species activity and the 
feedback loop of the HSR network (cf . Figure [2| . 

Figure [4] depicts the thermotolerance phenomena in the deterministic 



HSR model. Duration of the memory of the first temperature perturbation 



can be tracked by the activity of the HSP, as depicted in Figure S.4 




- 42. (on) 



3 



- 37. (off) 



300 400 
Time (min) 



Fig. 4: Thermotolerance in the heat-shock response: the substrate activity (solid) during 
the two consecutive heat-shocks (dotted) of 5°C over the homeostatis level of 37°C. The 
intensity of intoxication resulting from the amount of substrate (coloured area) depends 
on the time gap between heat-shocks. Interestingly, activity of the substrate in the second 
shock can be even higher than activity in the first shock, as shown for the time gap of 
240 min. 



In the stochastic model we introduce approximate perturbations as 



an independent, n-level Poisson process (cf. Figure S.3 see Section 5.1 



for details). This allows to stay within the same mathematical model 
and seamlessly perform stochastic simulations or model checking of the 
CTMC. 
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We define the notion of the thermotolerance during n-th heat-shock 
(n > 1) as the desensitisation coefficient: 

w 

where n-th response 7?. n is defined as: 

TZ n = max {#5 (t) - #5*} , (5) 

where #5* = Ejr (#<S) is a mean value in a steady state 7r, t n is a n-th 
heat-shock start time (we assume i n +i = oo if not specified otherwise), 
and the first response, by assumption, satisfies 1Z\ > 0. Such response 
measure represents the toxicity of the heat-shock: the higher the response 
the more likely the cell will die. For the deterministic model the species 



amount is simply a scaled value of ODE variable, corresponding to the 



mean value of a stochastic process random variable. Note also that the 
heat-shock response lZ n , as defined in Eq. (|5]), depends on the heat-shock 
duration At n ; if At n is short then the lZ n value may be equal to a value 
at the end of the heat-shock event (cf. Figure [4]). 

Figure [5] depicts value of the desensitisation coefficient T>2 for the sub- 
strate species, with respect to the time gap between heat-shocks. After the 
first heat-shock, at the time gap of the approximated memory loss, i.e. at 
time point t = t\ + At\ + 400 min, system is very close to the homeosta- 
sis steady state (data not shown). The slightly positive final level of T>2 m 
the stochastic model, as well as the overall difference with respect to the 
deterministic model may be attributed to the stochastic noise (we take 
maximum amount in Eq. (|5| 



3.2 Hyperthermia in multimodal oncological strategies 

Although hyperthermia exhibits a clear cytotoxic effect its efficacy 
is not enough to replace any one of the established therapy modalities 
when applied alone, but, undoubtedly, it is suitable enough to enhance 
the cell-killing effect of cytotoxic drugs or radiation [ "thermal chemosen- 
sitisation", "thermal radiosensitisation" [9]. In order to improve the ef- 
ficacy of anti-cancer therapies, it has been recently investigated how to 
combine different methods of cancer treatment into more effective multi- 
modal oncological strategies. Particular attention has been paid to treat- 
ments that involve hyperthermia as an adjuvant protocol for both radio 
and chemotherapy. A synergistic interaction between hyperthermia and 
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Fig. 5: The desensitisation coefficient T>2 for the substrate in the ODE model (black line) 
and in the |CTMC| plotted against the time gap between end of the first heat-shock and 
the beginning of the second heat-shock. Duration of both heat-shocks At n {n = 1,2) is 
equal to 71 min. Memory of the first heat-shock is lost when the desensitisation coefficient 
value stabilises around 0, which is approximately at 400 min for both mathematical models. 
Mean (yellow line) and standard deviation (orange line) of T>2 was calculated at selected 
time points (dots). Both estimators have a confidence interval with 95% confidence level. 
In case of the mean value the confidence interval width is less than 5- 10~ 3 , whilst for the 
stand ard devi ation the confidence interval is depicted as a strip. Estimators were calculated 
with 10 4 and 5 • 10 4 independent simulation samples for the first and the 



APMC 



using 

second moment respectively (see text for details). 



radiotherapy as well as various cytotoxic treatments has already been 
validated in pre-clinical studies [29 . Despite experimental evidences the 
precise mechanism of these interactions is not clear. 

It has been hypothesised that hyperthermia engages |HSR mechanism 
and its capacity, especially in cancer cells, is strongly limited, thus en- 
hancing the toxicity induced by a second modality [IT]. This synergistic 
effect of hyperthermia and other cancer therapies can be attributed to 
the much higher accumulation of denatured proteins (substrate), which 
are deadly for cell. To give a proof-of-concept modelling approach we 
investigate, by means of our deterministic HSR model, the temperature 
dependence of the heat-shock response in combination with bortezomib- 
induced inhibition of proteasome. 
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In order to incorporate into the model the inhibitory effect of borte- 



zomib we limit the HSP -assisted degradation of denaturated proteins 



(Eq. (r9)) and degradation of HSP itself (Eq. (|r8|)). More precisely, we 



linearly scale the reaction rate constants kg and kio, i.e. (1 — I) • k*, for 
I G [0,1], where I represents current inhibition level (when no drug is 
administrated 1 = whereas in case of maximum inhibition I = Iioo)- 

We used bortezomib pharmacodynamics as modelled by Sung and 
Simon [25] . Namely, the inhibition level linearly raises up to its maximum 
level at ii 100 =60 min, after which it decays with a half-life t\ 50 = 12 • ii 100 , 
i.e.: 

f -r^- for t < ilmn 



>) for t > tj 



1110 



where ki = ln(2)/ti 50 « 9.63- 10" 4 . The maximum inhibition level Iioo di- 
rectly corresponds to the drug dose. For a maximum tolerated bortezomib 
dose, Iioo is equal to ca. 65%, while for some of the next-generation pro- 
teasome inhibitors, such as carfilzomib or ONX-0912 both of which are 
in clinical development, it was possible to reach over 80% proteasome 
inhibition in blood [with consecutive-day dosing schedules: 116] . 

Figure [6] depicts substrate activity, which defines level of heat-shock 
response TZ% (cf . Eq. ([5| ) , with respect to a unimodal moderate hyperther- 
mia treatment and a full range of unimodal proteasome inhibition treat- 
ment, as represented by substrate activity. Please note that 1Z± is limited 



by the saturation limit of S, which is bounded by P to t (cf- Eq. (cl )). 

When both therapies are simulated simultaneously, substrate accumu- 
lation, indeed, reaches higher level than it would be possible with applica- 
tion of only one of the treatment modalities. Figure [7] depicts this syner- 
gistic effect. We observed monotonic increase in response with respect to 
both temperature and maximum inhibition level for all heat-shock appli- 
cation times up to ti 100 (not shown). We found that multimodal toxicity 
response can increases by over 40% with respect to a unimodal hyperther- 
mia response for a maximum inhibition level equal to reported 65%, up 
to over 80% increase for a theoretical maximum of 100% of proteasome 
inhibition. Moreover, we found that optimal time to start hyperthermia 
treatment in combination with bortezomib application is t* ~ 38 min. 
Interestingly, this is not in agreement with a maximum area under the 
bortezomib inhibition curve (AUC), a common pharmacokinetic efficacy 
measure. For the heat-shock duration At\ = 71 min AUC is reached at 



t\ ~ 56 min (see Figure S.5). Timing of heat-shock in the optimal mul- 



timodal treatment strategy t\ can be explained by two observations (cf. 
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Figure [6]) . Firstly, time required for denaturated proteins to peak after 
the beginning of a heat-shock is roughly the same as the time gap be- 
tween t* and tj 100 (22 min). Secondly, at t\ the inhibition itself has still 
a relatively low effect. This way, inhibition peak coincides with the period 
of maximum temperature-induced toxicity, at which HSR mechanism is 
the most occupied, thus, effecting in the optimal synergistic toxicity. 



4 Conclusions and Discussion 



We demonstrated feasibility and practical potential of the PMC tech- 
nique, more specifically its lesser known approximate variant, in a case 
study of analysis of the heat-shock response mechanism. We compara- 
tively studied its deterministic and stochastic variants, including the in- 
vestigation of the thermotolerance phenomenon. We found deterministic 
approach to be a valid approximation of stochastic model, i.e. both vari- 
ants presented very similar results with respect to the stochastic mean. 
However, the stochastic model presented a constant high standard devi- 
ation of the thermotolerance phenomenon (ca. 20% of its expected max- 
imum level). 



Next, by mechanistic modelling of HSR we were able to support the 
common belief that the combined cancer treatment strategies can more 
effectively increase cytotoxicity of denatured proteins in cancer cells than 
unimodal strategies. Moreover, we presented an optimal starting time for 
a moderate hyperthermia treatment in combination with a proteasome 
inhibitor application. This is an example of how mechanistic modelling 
can enrich pharmakokinetic measures of optimal efficacy, such as AUC 
which basically is an optimisation only with respect to system's input. 

We suggest that the synergistic effect of hyperthermia and other can- 
cer treatment modalities (like chemo- and radiotherapies) is caused by in- 
creased accumulation of denatured proteins, i.e., heat and drug-sensitive 
proteins or heat and radiation sensitive proteins. This results in bigger 
demand for heat-shock proteins and higher selective barrier for cells. 

Our model-based analysis proves successful in reproducing experimen- 
tal knowledge of key aspects of hyperthermia treatment, and as such of- 
fers reasonable framework for studying its connections with heat-shock 
response. However, all of the kinetic models of molecular biological sys- 
tems and means of their analysis are incomplete due to constraints under 
which these models are formulated. In this regard, we would like to point 
out that this is a proof of concept model-based analysis and there are 
many issues to address within this work. For instance, we omitted the 
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investigation of the day-based strategies of multimodal treatment. This 
is because we found that in our model the single cell level thermotoler- 
ance duration (ca. 6.5h) is much shorter than the bortezomib decay (12h 
half-life), thus, making the latter a determining factor for a standard, 
daily dosing schedule. The inconsistency between reported (24-48h) and 
simulated thermotolerance duration is one of the open modelling issues. 
With regard to a function determining the protein denaturation rate con- 
stant (Eq. 0), we had to adjust its constant not only to time scale but 
also to fit the original, experimental data (as pointed out by Andrzej 
Mizera in personal communication). Last but not least, to provide solid, 
quantitative results, our model requires more extensive calibration with 



respect to experimental data (cf. Figure S.2), including the temperature- 
dependent saturation behaviour [cf. [26] . Nevertheless, undoubtedly, our 
proof of concept analysis presents a valuable advance towards model- 
based understanding of hyperthermia treatment strategies. 

5 Methods 

5.1 Probabilistic model checking 



PMC requires two inputs: a description of the probabilistic system, 
usually given in some high-level modelling language; and a specification of 
desired properties of the system, typically in probabilistic temporal logics 



such as: probabilistic computation tree logic PCTL| 7 j for probabilistic 
timed automata or IPCTLf for discrete time Markov chain and Markov 
decision process, both extensions of CTL*; or continuous stochastic logic 



|CSL| 2], an extension of |PCTL| for |CTMC| Adopting probabilistic model 



checking paradigm allows us to take advantage of a efficient probabilistic 
model checker PRISM, which implements all of the mentioned |PMC| 
variants. 

Models are described using the compositional PRISM language. It is 
a simple, state-based language, supporting a wide range of probabilistic 
models and their extensions with costs and rewards. In our experiments 



we have used CTMC as our principal model. 

In PRISM, properties are expressed in temporal logic suitable for 
chosen type of the probabilistic model. PRISM incorporates extensions 
to these logics for quantitative specifications and rewards [see, e.g., [TO] . 
In our experiments we have used |CSL| for |CTMC model checking, An 
example of a typical property is: 



P=? (0) 



(7) 
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It asks about the probability that a path in the model state space sat- 
isfies (p. In other words, PRISM computes the fraction of all paths over 
model states that satisfy <f>. 



Approximate PMC (APMC) 

The model analysed in our experiments is so large that standard model 
checking becomes unfeasible. We have thus decided to use more efficient 
techniques that trade accuracy for efficiency. PRISM includes a range of 
approximate analysis techniques, based on sampling: generating a number 
of random paths through the model state space (model trajectories) by 
stochastic simulation, and evaluating the property on each run [for an 
overview see 18 . This information is used to compute an approximate 
result. In our experiments we have used the simplest and the most flexible 
CI implementation of APMC [cf. 6 . CI computes a confidence interval 
for the value of a given property, such as the one in Eq. @. CI method 
takes three parameters: interval half-width w, confidence level a and the 
number of sampled paths N. On the basis of evaluation of sampled paths, 
the method computes an approximate value y of a formula, determining 
the confidence interval (y — w,y + w). The exact (unknown) value x of 
the formula falls into the confidence interval with probability 1 — a. 



5.2 Estimators 

Mean value of amount of each species was estimated using the confi- 



dence interval APMC method to verify the reward-based property: 

R{#s = ?} (I = burn-in time) , 
where #5 reward for each species S is defined as: 
rewards "#S" true : S; endrewards 

We start the stochastic process with a single point distribution, according 
to RRE values (cf. Figure [3]), thus the burn-in time, which we assumed 



300 min for T = 37°C and 600 min for T = 42°C. 

It is impossible to query for higher central moments in PRISM in 
a single run. Therefore, the variance value of amount of each species was 
estimated from the unbiased mean value and second moment estimators, 
i.e. 
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Second moment E (#S 2 ) of species amount was estimated analogously to 
the mean value, using the confidence interval APMC method (see above). 
Having symmetric confidence intervals: 



E(X) € (E(X)-ai,E(X) + aiJ and 
E (X) € (eJx*) - a 2 , e7^2) + a 2 ) , 



with a confidence level, the unbiased moments-based variance estimator 
Var (X) has an asymmetric confidence interval: 



Var (X) := (e[X^) - a 2 ) - (e(X) + a x 

< Var"(X) < 



2 



E(X 2 ) + a 2 ) - (E(X) -at) =:Var(X). 

(8) 



Analogously, for an unbiased VMR estimator 
VMR (X) = Var (X) ^E (X) , we get asymmetric confidence interval, with 
a confidence level, from the following inequalities: 

Var (X) Var ( X) 



VMR(V) := < VMR(V) < v ' =: VMR(X). 

E (X) +ai E (X) - ax 

For a symmetric confidence interval with a confidence level (cf. Table [6]), 
we simply take: 

VMlTpO t max (vMRpf ) - VMR (V),VMR (X) - VMRpT)) . 



The T>2 mean value for heat-shocks time gap t\ — (ii + At\) is esti- 
mated using the confidence interval APMC method to verify the reward- 
based property: 

:= R{^ =? } (I = 1-05 • (t 2 + At 2 )) 

where for i = 1,2, first and second moment T> 2 rewards are defined as: 
rewards 

"V\" true: {S^ - S*)/(S^ ax - S*) ; 
n V\" true: ( (Sf ax - S*)/(Sf ax - S*) f ; 
endrewards 
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Here, Sf ax is an additional stochastic model variable, which is a witness of 
the maximum value of the substrate variable S, during and after the z-th 
heat-shock. Introduction of such variable in PRISM modelling language 
can be d one seamlessly, i.e. without affecting the behaviour of the original 
Finally, we have E(S) = 1 - as E (1 - X) = 1 - E (X), 



CTMC 



and Var(£> 2 ) = (f> 2 - {(f) 1 ) 2 , as Var(l-X) = Var(X). The unbiased 




standard deviation estimator SD (X^) = V Var (2^), nas the following a 
level confidence interval: 



Var (X) < SD (X) < yj Var (X), 

provided that the variance estimator's precision is high enough, i.e. 
Var(X) > (cf. Eq. d8b). 



5.3 Approximate stochastic perturbation strategy 

In our experiments we model scenarios of consecutive heat-shocks im- 
posed on the model, separated with gaps. In principle, we control time 
points when the heat shocks are activated or inactivated. However, to 
do that in a probabilistic model one has to either save the whole distri- 
bution over a possibly infinite state space or, if analysis is based solely 
on stochastic simulations, modify the simulation algorithm to test again 
a current time point and queued deterministic events. Essentially, these 
are solutions based on moving to a more general type of models, like 
Markov decision processes, but the price to pay would be complication 



and impracticability of analysis (e.g., the APMC does not handle Markov 
decision processes). 

We have thus chosen to consider heat-shock as an approximate random 
perturbation event and, thereby to stay within the same mathematical 
model and seamlessly perform stochastic simulations or model checking 
of lCTMCI 

Our approach relies on introduction of the sequence of rij-counting 



Poisson processes (a special case of one-dimensional CTMC), indepen- 
dent of other state variables. Each time i-th process reaches value rtj, it is 
replaced by another nj+i-counting Poisson processes with rate m+i -Ti+i, 
where l/r^+i = ij+i — (U + Ati), i.e. r^+i is an inverse of a time gap be- 
tween perturbations i and Because time between consecutive Poisson 
process events is exponentially distributed, i.e. Tij ~ Exp (n,rj), the ex- 
pected time for an approximate perturbation event equals to the time 
gap between deterministic perturbations time, i.e. E (X^=i^ij) = l/ r «- 
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Moreover, due to independence of a time of occurrence of each count of 
the Poisson process Tij, variance of the i-th perturbation event, calculated 

independently of prior i — 1 perturbations, is equal to Var ^X^j=i ^Vj = 

Sj=i Var (T^) = ijuirf . In other words the precision of single pertur- 
bation event is linearly proportional to the number of counting levels and 
inverse linearly proportional to the time of its occurrence. 

In PRISM modelling language we introduce two heat-shock events, 
i.e. four perturbation events of temperature parameter T, using Poisson 
processes with a common n levels for each parameter and an additional 
perturbation number counter i. Using the compositional description of 
variables, which represent CTMC| state, and commands, which change 



CTMC state, the independent perturbation events module (i.e. not syn- 



chronised with any other commands) is defined as: 



ctmc 



const double tl; // time offset for 1st heat-shock 
const double tdl; // duration of 1st heat-shock 
const double t2; // time offset for second heat-shock 
const double td2; // duration of 2nd heat-shock 
const int Td; // heat-shock temperature delta 
const int n; // event switcher levels 



module events 

i : [0. .4] init 0; // number of perturbation 

ps: [0..1] init 0; // perturbation switcher 

cnt : [l..n] init 1; // actual poisson process variable 

//pre 1st heat-shock 

] i = & cnt < n -> n/tl: (cnt ,= cnt+l) ; 

] i = & cnt = n -> n/tl: (i'=l)&(ps'=l)&(cnt'=l) ; 

//1st heat-shock 

] i = 1 & cnt < n -> n/tdl: (cnt ' =cnt+l) ; 

] i = 1 & cnt = n -> n/tdl: (i , =2)&(ps , =0)&(cnt , =l) ; 

//pre 2nd heat-shock 

] i = 2 & cnt < n -> n/t2: (cnt ,= cnt+l) ; 

] i = 2 & cnt = n -> n/t2: (i , =3)&(ps , =l)&(cnt ) =l) ; 
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//2nd heat-shock 

[] i = 3 & cnt < n -> n/td2: (cnt'=cnt+l) ; 

[] i = 3 & cnt = n -> n/td2: (i , =4)&(ps , =0)&(cnt , =l) ; 
endmodule 



formula T = (37+Td*ps) ; // temperature for misfolding rate 

Clearly, increasing n increases the number of states of the model, and 
thus makes the experiments less efficient. A suitable value of the parame- 
ter n has been chosen experimentally by visual assessment of the substrate 
stochastic trajectories precision, taking under consideration simulations 



efficiency; see Figure S.3 caption for details on choice of n in our HSR 
model. 



5.4 Additional software tools 



We defined model using the SBML-shorthand notation. The |RRE" 



model has been numerically solved using the MathSBML package 
of Mathematica software. The latter was used for creation of most of 
the plots and for collection of results. To generate the PRISM model, we 
have used a prototype SBML translator which generates model specifica- 
tion in the PRISM language. Some minor adjustments, such as factori- 



sation of parameters or accounting for mass conservation laws (cl)-(c3), 
have been done manually. The approximate stochastic perturbation strat- 
egy has also been encoded manually in the PRISM model according to 
presented scheme. All stochastic simulations and the confidence interval- 
based [SPMU] were done using PRISM. 
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Supplementary Material 
S.l Model file 



HSR| model analysed in this paper is available as a XML file in the SBML 



format (deterministic version), as well as a text file in the PRISM model 
format (stochastic version). 



S.2 Figure: HSE:HSF3 fit to the experimental data 




0.0 - 

50 100 150 200 250 300 350 

Time (min) 



Fig. S.2: The HSE:HSF 3 fit to the experimental data [TJ Fig. 8A], for a constitutive 42°C 
heat-shock. We assumed relative 15% error plus 2% of the peak value. The resulting fit 
gives x 2 (10)=2.49 (with p-value equal to 99%), for concentration scale equal to 7.3 x 10 -3 
"levels" to a.s.M. 



S.3 Figure: Approximate stochastic perturbation of 
temperature in the HSR model 
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Fig. S.3: Comparison of simulations of approximate perturbations in the stochastic 
model with respect to ODE simulations. The comparison is shown for different levels n 
of the independent Poisson process, which approximates the deterministic scheme of 
perturbations of the temperature (left column). Plots show examples of 50 stochastic 
simulations, and mean +/— standard deviation, estimated from 1000 samples, for the 
temperature perturbation scheme (left column) and the substrate amount (right column). 
Variance of time of a temperature perturbation event T is a function of an perturbations 
time gap t and the precision levels n, i.e. Var(T) = t 2 /n. Therefore, the mam time 
approximation error for our perturbation scheme lies within an event of the second 
heat-shock with t = 343 min (cf. left column). Note though, that errors for subsequent 
perturbation times accumulate, so the biggest error can be observed for the end of the 
second heat shock. Nevertheless, because start of the second heat-shock creates high 
amounts of molecules in the system, it is the most error influencing event with respect 
to species amounts expected from the |RRE| numerical solutions (cf. right column). We 
found that for n = 2 13 the approximation of mean and standard deviation of the substrate 
amount in the second heat-shock is in a good agreement with much more precise values 
in the first heat-shock and with IRREl solutions. 



Figure on the next page. 
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S.4 Figure: The substrate and HSP heat-shock response 
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S.5 Figure: area under inhibition curve with respect to 
heat-shock time 
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Fig. S.5: Area under the bortezomib inhibition curve (AUC) versus heat-shock time ti, 
for a fixed heat-shock duration Al\ — 71 min. More specifically, AUC( tl . tl+Z } tl ) = 
Iti +Atl ' M/'ioo dt; maximum value is reached for t\ — 56.24 min. 



Bibliography 



[1] Abravaya, K., Phillips, B., Morimoto, R.I.: Attenuation of the heat 
shock response in HeLa cells is mediated by the release of bound heat 
shock transcription factor and is modulated by changes in growth and 
in heat shock temperatures. Genes & Development 5(11), 2117-2127 
(1991) 

[2] Aziz, A., Sanwal, K., Singhal, V., Brayton, R.: Verifying continuous 
time markov chains. In: Proceedings of the 8th International Con- 
ference on Computer Aided Verification (CAV'96). LNCS, vol. 1102, 
pp. 269-276. Springer (1996) 

[3] Calder, M., Gilmore, S., Hillston, J., Vyshemirsky, V.: Formal meth- 
ods for biochemical signalling pathways. In: Boca, P., Bowen, J. P., 
Siddiqi, J. (eds.) Formal Methods: State of the Art and New Direc- 
tions, pp. 185-215. Springer London (2010) 

[4] Calder, M., Vyshemirsky, V., Gilbert, D., Orton, R.: Analysis of sig- 
nalling pathways using continuous time Markov chains. In: Priami, 
C, Plotkin, G. (eds.) Transactions on Computational Systems Bi- 
ology VI, LNCS, vol. 4220, pp. 44-67. Springer Berlin / Heidelberg 
(2006) 

[5] Charzyhska, A., Nalecz, A., Rybihski, M., Gambin, A.: Sensitivity 
analysis of mathematical models of signaling pathways (2012), in 
press (BioTechnologia) 

[6] Grosu, R., Smolka, S.A.: Monte Carlo model checking. In: Proceed- 
ings of the 11th International Conference on Tools and Algorithms 
for the Construction and Analysis of Systems (TACAS'05). LNCS, 
vol. 3440, pp. 271-286. Springer (2005) 

[7] Hansson, H., Jonsson, B.: A logic for reasoning about time and reli- 
ability. Formal Aspects of Computing 6(5), 512-535 (1994) 

[8] Heath, J., Kwiatkowska, M., Norman, C, Parker, D., Tymchyshyn, 
O.: Probabilistic model checking of complex biological pathways. 
Theoretical Computer Science 391(3), 239-257 (2008) 

[9] Hildebrandt, B., Wust, P., Ahlers, O., Dieing, A., Sreenivasa, G., 
Kerner, T., Felix, R., Riess, H.: The cellular and molecular basis 
of hyperthermia. Critical Reviews in Oncology /Hematology 43(1), 
33-56 (2002) 

[10] Kwiatkowska, M., Norman, G., Pacheco, A.: Model checking ex- 
pected time and expected reward formulae with random time 



34 

bounds. Computers & Mathematics with Applications 51(2), 305 — 
-316 (2006) 

[11] Kwiatkowska, M., Norman, G., Parker, D.: Using probabilistic model 
checking in systems biology. ACM SIGMETRICS Performance Eval- 
uation Review 35(4), 14-21 (2008) 

[12] Lepock, J.R., Frey, H.E., Ritchie, K.P.: Protein denaturation in in- 
tact hepatocytes and isolated cellular organelles during heat shock. 
The Journal of Cell Biology 122(6), 1267-1276 (1993) 

[13] Michaelis, L., Menten, M.: Die Kinetik der Invertinwirkung. Bio- 
chemische Zeitschrift 49, 333-369 (1913) 

[14] Milo, R., Jorgensen, P., Moran, U., Weber, C, Springer, M.: BioN- 
umbers — the database of key numbers in molecular and cell biology. 
Nucleic Acids Research 38(Database issue), D750-753 (2010) 

[15] Mizera, A., Gambin, B.: Stochastic modelling of the eukaryotic heat 
shock response. Journal of Theoretical Biology 265(3), 455-466 
(2010) 

[16] Molineaux, S.M.: Molecular pathways: targeting proteasomal protein 
degradation in cancer. Clinical cancer research: an official journal of 
the American Association for Cancer Research 18(1), 15-20 (2012) 

[17] Neznanov, N., Komarov, A. P., Neznanova, L., Stanhope-Baker, P., 
Gudkov, A. V.: Proteotoxic stress targeted therapy (PSTT): induc- 
tion of protein misfolding enhances the antitumor effect of the pro- 
teasome inhibitor bortezomib. Oncotarget 2(3), 209-221 (2011) 

[18] Nimal, V.: Statistical Approaches for Probabilistic Model Checking. 
Master's thesis, Oxford University Computing Laboratory (2010) 

[19] Owens, N., Timmis, J., Greensted, A., Tyrrell, A.: Modelling the 
tunability of early t cell signalling events. In: Bentley, P., Lee, D., 
Jung, S. (eds.) Artificial Immune Systems, LNCS, vol. 5132, pp. 12- 
23. Springer Berlin / Heidelberg (2008) 

[20] Peper, A., Grimbergen, C.A., Spaan, J. A., Souren, J.E., van Wijk, 
R.: A mathematical model of the hsp70 regulation in the cell. Inter- 
national Journal of Hyperthermia: The Official Journal of European 
Society for Hyperthermic Oncology, North American Hyperthermia 
Group 14(1), 97-124 (1998) 

[21] Petre, I., Mizera, A., Hyder, C.L., Meinander, A., Mikhailov, A., 
Morimoto, R.I., Sistonen, L., Eriksson, J.E., Back, R.: A simple 
mass-action model for the eukaryotic heat shock response and its 
mathematical validation. Natural Computing 10(1), 595-612 (2010) 

[22] Proctor, C.J., Lorimer, I. A. J.: Modelling the role of the Hsp70/Hsp90 
system in the maintenance of protein homeostasis. PLoS ONE 6, 
e22038 (2011) 



35 



[23] Rieger, T.R., Morimoto, R.L, Hatzimanikatis, V.: Mathematical 
modeling of the eukaryotic heat-shock response: dynamics of the 
hsp70 promoter. Biophysical Journal 88(3), 1646-1658 (2005) 

[24] Shapiro, B.E., Hucka, M., Finney, A., Doyle, J.: MathSBML: a pack- 
age for manipulating SBML-based biological models. Bioinformatics 
(Oxford, England) 20(16), 2829-2831 (2004) 

[25] Sung, M.H., Simon, R.: In silico simulation of inhibitor drug effects 
on nuclear factor- kappaB pathway dynamics. Molecular pharmacol- 
ogy 66(1), 70-75 (2004) 

[26] Szymahska, Z., Zylicz, M.: Mathematical modeling of heat shock pro- 
tein synthesis in response to temperature change. Journal of Theo- 
retical Biology 259(3), 562-569 (2009) 

[27] Wegele, H., Miiller, L., Buchner, J.: Hsp70 and hsp90~a relay team 
for protein folding. Reviews of Physiology, Biochemistry and Phar- 
macology 151, 1-44 (2004) 

[28] Wilkinson, D.J.: Stochastic Modelling for Systems Biology. CRC 
Press (2011) 

[29] Wust, P., Hildebrandt, B., Sreenivasa, G., Rau, B., Gellermann, J., 
Riess, H., Felix, R., Schlag, P.M.: Hyperthermia in combined treat- 
ment of cancer. The Lancet Oncology 3(8), 487-497 (2002) 

[30] van der Zee, J.: Heating the patient: a promising approach? An- 
nals of oncology: official journal of the European Society for Medical 
Oncology / ESMO 13(8), 1173-1184 (2002) 



